% This code simulates the model for EIS = 1 under RE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SIMDATA = Sim_EIS_1_RE(l_q , l_q_burn , N_sim , c_0 , Var_Names)

% INPUT:
% -- l_q: Number of sample periods
% -- l_q_burn: Number of burn-in periods
% -- N_sim: Number of simulations
% -- c_0: starting value of log real per capita consumption level
% -- Var_Names: names of variables to be stored
%
% OUTPUT:
% -- SIMDATA: stored simulations as a cell

%% 1. LOAD PARAMETERS
I = param_EIS_1(0);

%% 2. SIMULATIONS
% Additional parameters for simulation
T = l_q;                %Number of sample periods
T_burn = l_q_burn;            %Number of burn-in periods
T_total = T_burn + T;     %Total periods including burn-in periods for simulation
d_0 = c_0 + I.d_c_mean;   %Starting value of log dividend

%% 2.1 Simulate the economy
% Simulate quarterly iid consumption growth, column for each simulation
cg_full = normrnd(I.mu , I.sigma , T_total , N_sim); 

% Simulate quarterly shocks to dividend growth
dg_shock = normrnd(0 , I.sigma_d , T_total , N_sim); 

% Log real consumption level from 1 to T_total
c_full = c_0 + cumsum(cg_full); 

% Log real dividend level from 1 to T_total
d_full = zeros(T_total , N_sim); 

d_full(1 , :) = I.lambda * cg_full(1 , :) + (1 - I.alpha) * d_0 + I.alpha * c_0 + I.alpha * I.mu_dc + dg_shock(1 , :);

for k = 2 : T_total
  d_full(k , :) = I.lambda * cg_full(k , :) + (1 - I.alpha) * d_full(k - 1 , :) + I.alpha * c_full(k - 1 , :) + I.alpha * I.mu_dc + dg_shock(k , :);
end

% Log dividend growth from 1 to T_total
dg_full = [d_full(1 , :) - d_0 ; d_full(2 : end , :) - d_full(1 : end - 1 , :)]; 

% Four-quarter change
cg_yoy_full = (c_full(5 : end , :) - c_full(1 : end - 4 , :))/4;
dg_yoy_full = (d_full(5 : end , :) - d_full(1 : end - 4 , :))/4;

% Construct State Variable \tilde{\mu}
tmu_full = I.mu * ones(T_total , N_sim);

% Calculate Risk-free rate from 1 to T_total
rf_full = -I.tmu_m + tmu_full - 0.5 * I.gamma^2 * I.sigma^2; 

%% 2.2 Calculate price
P_full = zeros(T_total , N_sim);                 %Price level

% Max number of terms to calculate explicitly
J = round(log(0.000001)/log(1 - I.alpha)); 

% Calculate the first J terms of A_{c,n}: A_{c,1},...,A_{c,J} as in Equation (D.17)
A_n_inc = ( (I.lambda - 1) * (1 - I.alpha).^(0 : J - 1) + 1 - I.gamma ).^2;
      
A_N = cumsum(A_n_inc); % These are the A_{c,n}

% Further incremental terms to A_{c,n} will be approximated by A_n_inc_lim
A_n_inc_lim = ( 1 - I.gamma )^2;

% Calculate the first J dividend strips' prices at time t. first calculate
% the multiplicative term that is invariant to time on the RHS of Equation (D.16)
p_cons = ( I.tmu_m * (1 : J) ... % n\tmu_m
       + (1 - (1 - I.alpha).^(1 : J)) * (I.mu_dc + (I.lambda - 1) * I.mu/I.alpha) ...%(1 - (1-\alpha)^n)\mu_{dc}
       + 0.5 * A_N * I.sigma^2 ... %0.5 * A_{c,n} \sigma^2
       + 0.5 * (1 - (1 - I.alpha).^(2 * (1 : J))) * I.sigma_d^2/(1 - (1 - I.alpha)^2) )'; % 0.5 * A_{d,n} \sigma_d^2

% Multiplier that is invariant with respect to simulations or time in
% Equation (D.31) and (D.32)
approx_cons = exp( I.mu_dc + (I.lambda - 1) * I.mu/I.alpha + 0.5 * I.sigma_d^2/(1 - (1 - I.alpha)^2) ...
            + (J + 1) * I.tmu_m + 0.5 * A_N(J) * I.sigma^2 + 0.5 * A_n_inc_lim * I.sigma^2 ) ...
            /(1 - I.delta);

% Calculate total price using Equation (D.31)        
for n = 1 : T_total 
    %Calculate p_t^n for maturity up to J
    p_t_n_temp = d_full(n , :) ... %d_t from LHS
  + bsxfun(@times , ( 1 - (1 - I.alpha).^(1 : J) )' , c_full(n , :) - d_full(n , :)) ...
  + p_cons;

    %Calculate aggregate price
    P_full(n,:) = sum(exp(p_t_n_temp) , 1) ... % First J dividend strip
                + exp(c_full(n , :)) * approx_cons;    
end

%Calculate trailing 4-quarter sum of dividends
D_ttm_full = movsum(exp(d_full) , 4 , 'Endpoints' , 'discard');

%% 2.3 Calculate other quantities

% P/D ratio
PD_ttm_full = P_full(4 : end , :) ./ D_ttm_full; 
pd_ttm_full = log(PD_ttm_full);

PD_full = P_full./exp(d_full);
pd_full = log(PD_full);

% Objective realized returns
ret_full = log( P_full(2 : end , :) + exp(d_full(2 : end , :)) ) ... % log(P_{t+1} + D_{t+1})
         - log( P_full(1 : end - 1 , :) ); % - log(P_t)

EXRET_full = exp(ret_full) - exp(rf_full(1 : end - 1 , :)); % R_{t+1} - R_{f,t}

exret_full = ret_full - rf_full(1:end-1,:); %r_{t+1} - r_{f,t}

%% 3. Simulated Data

%% 3.1 Drop burn-in periods
% All data are re-labled from 1 to T as realizations. 
% Be careful with risk-free rate
%-------------------------------------
% Endowment Process
%-------------------------------------
cg = cg_full(end - T + 1 : end , :);
dg = dg_full(end - T + 1 : end , :);
c  =  c_full(end - T + 1 : end , :);
d  =  d_full(end - T + 1 : end , :);
cg_yoy = cg_yoy_full(end - T + 1 : end , :);
dg_yoy = dg_yoy_full(end - T + 1 : end , :);

%-------------------------------------
% State variable proxies
%-------------------------------------
%tmu  =  tmu_full(end - T + 1 : end , :);
tmu = nan(T , N_sim);
tmud = nan(T , N_sim);
tmur = nan(T , N_sim);
LTG = nan(T , N_sim); 
Y1G = nan(T , N_sim);

%-------------------------------------
% Price Quantity
%-------------------------------------
rf = rf_full(end - T + 1 : end , :);

P = P_full(end - T + 1 : end , :);
D_ttm = D_ttm_full(end - T + 1 : end , :);
pd_ttm = pd_ttm_full(end - T + 1 : end , :);
pd = pd_full(end - T + 1 : end , :);

ret = ret_full(end - T + 1 : end , :);
EXRET = EXRET_full(end - T + 1 : end , :);
exret = exret_full(end - T + 1 : end , :);


%% 3.2 Store simulations
SIMDATA = cell(length(Var_Names) , 1);
for i = 1 : length(Var_Names)
    SIMDATA{i} = eval(Var_Names{i});
end

end